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Abstract:  The  propagation  of  shear  (S)  and 
compression  (P)  waves  within  the  earth  allows 
geologists  to  track  seismic  events  and  to  identify 
subterranean  stmcture.  Highly  specialized 
geological  based  computer  programs  developed 
have  been  instmmental  in  determining  the 
location  and  characteristics  of  natural 
phenomena  (e.g.,  earthquakes)  and  man-made 
activity  (e.g.,  nuclear-blast  tests).  Use  of  these 
internally  developed  programs  requires  regular 
maintenance,  and  the  reliance  on  supercomputers 
limits  broad  accessibility.  This  paper  seeks  to 
demonstrate  that  commercially  available 
software  mnning  on  desktop  computational 
resources  can  provide  accurate  solutions  to  an 
important  subset  of  problems  associated  with 
wave  propagation  in  the  geophysical  domain. 

The  work  presented  here  uses  COMSOL 
Multiphysics  to  solve  the  equilibrium  equations 
for  a  time-varying  system  using  the  finite 
element  method.  This  work  focuses  on 
developing  a  benchmark  solution  of  a 
homogeneous  half-space  loading  with  an  impact 
and  develops  a  general  closed- form  solution 
against  which  to  compare  the  computational 
results.  These  results  show  the  ability  to  resolve 
both  S  and  P  wave  across  the  computational 
domain.  Thus,  COMSOL  Multiphysics  running 
on  desktop  computational  resources  provides 
sufficiently  accurate  results  for  critical 
geophysical  wave  propagation  problems. 
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1.  Introduction 

Scientists  and  engineers  that  seek  to  understand 
elastic  wave  propagation  in  geological  structures 
typically  consider  the  generic  problem  of  a 
seismic  wave  that  is  generated  at  a  source, 
propagates  through  a  media,  and  is  measured  at  a 
receiver.  Some  researchers  are  concerned  with 
natural  sources,  e.g.  earthquakes,  while  others 
focus  on  man-made  sources,  e.g.  explosions. 


Some  who  study  these  areas  seek  to  understand 
the  nature  of  the  source,  while  others  use  well 
characterized  sources  and  seek  to  characterize 
the  medium.  These  researchers  may  be  seeking 
to  predict  future  earthquakes,  locate  natural 
resources,  or  identify  and  locate  specific  human 
activities.  However,  in  all  these  cases,  the 
overarching  physics  of  elastic  wave  propagation 
in  a  solid  medium  remain  the  same. 

Seismic  waves  propagating  in  bulk  material  do 
so  as  either  compressional  waves,  P,  where 
material  translates  in  the  direction  of  wave 
propagation,  and  shear  waves,  S,  where  material 
translates  perpendicular  to  the  direction  of  wave 
propagation.  P- waves  travel  faster  than  S-waves. 
Another  important  class  of  waves  is  surface 
waves;  these  waves  develop  due  to  an  energy 
concentration  near  the  Earth’s  surface  and 
consequently  propagate  in  two  dimensions.  As 
such,  surface  waves  decay  as  1/r  while  body 
waves  decay  as  1/r2  due  to  their  propagation  in 
three-dimensions.  Thus,  energy  measured  at  a 
near-surface  receiver  located  a  significant 
distance  from  the  source  is  dominated  by  surface 
waves.  Typically,  these  types  of  waves  are 
referred  to  as  either  Love  or  Rayleigh  waves. 
Love  waves  are  shear  waves  that  have  been 
polarized  in  the  horizontal  direction  (parallel  to 
the  Earth’s  surface)  while  Rayleigh  waves  are  a 
mixture  of  P  waves  and  S  waves  that  have  been 
polarized  in  the  vertical  direction. 

Geophysics  studies  typically  use  specialized 
computer  software  running  on  supercomputers  to 
simulate  wave  propagation  in  geophysical 
domains.  This  work  demonstrates  that  a  finite 
element  based  software  package  that  is 
commercially  available,  COMSOL  Multiphysics, 
solves  these  wave  propagation  problems  using 
readily  available  desktop  computational 
hardware. 
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2.  Method 

This  work  uses  the  Structural  Mechanics  Module 
available  in  COMSOL  Multiphysics  to  develop 
computational  models  for  a  range  of  problems. 

These  models  increase  with  complexity  ranging 
from  a  simple  point  source  in  an  infinite  solid  to 
a  volume  source  that  represents  experimentally 
measure  forces  applied  in  a  layered  half-space. 

2.1  Point  Source  in  an  Infinite  Media 

Initially,  this  work  focuses  on  modeling  a  point 
source  located  in  an  infinite  domain.  To  provide 
a  point  of  comparison  for  the  finite  element 
models,  an  analytical  solution  for  this  problem 
was  developed,  and  is  described  as  follows. 
Using  the  notation  Uy  to  represent  the 

displacement  in  the  i  -  direction  due  to  a 
concentrated  point  force,  specified  as 
f  (x,  t;  %)  =  f0  (t)s(\  -fyj,  applied  in  the 

j  —  direction  at  the  point  £, ,  the  displacements 

over  the  domain  may  be  determined  by  solving 
the  following  equation; 


This  work  represents  the  typical  source  force  as: 
/0(?)  =  H{t)  =  crte  ,  where  H(t  )  is  the 
Heaviside  step  function,  a  is  a  parameter 
controlling  duration  and  amplitude  of  the  source, 
for  additional  information  see  [1]. 

Three  finite  element  models  of  a  spherical 
domain  loaded  by  a  point  source  were  developed 
in  COMSOL  Multiphysics.  The  volume  was 
represented  in  three-dimensions  using  a  quarter 
symmetric  finite  element  model.  To  reduce  the 
computational  cost  of  this  model,  two- 
dimensions  models  where  constructed  using  an 
axisymmetric  and  a  plane  strain  formulation. 

2.2  Point  Source  in  a  Semi-Infinite  Media 

To  extend  this  work  to  include  surface  wave,  a 
solution  to  Lamb’s  problem  of  a  source  loading 
at  the  surface  of  a  semi-infinite  domain  is 
developed.  Again,  a  closed  form  solution  was 
developed  as  a  point  of  comparison  for  the  finite 
element  solution  of  this  problem.  In  this  closed 
form  solution,  the  vertical  stress  acting  on  the 
free  surface  is  denoted  as  ; 

where 


4  —  frf0(t-r)dr 
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fit)  =  m2cr33(t)  denotes  the  total  applied 
force 

and 


where 

p  is  the  material  density, 
r  —  |x  —  is  the  distance  from  the  source 
x .  -  £ 

Yi  — -  are  the  direction  cosines  of 
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vector,  and 
a  =  cp 
and 
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V  P 


P  =  CS 


represent  the  speeds  of  pressure  and  shear 
waves,  respectively.  In  these  equations, 

/L,  p  are  Lame  elastic  constants, 


a  denotes  the  radius  of  the  circle  over  which 
the  stress  is  applied. 

From  this  loading,  the  vertical  displacement  is 
denoted  w(t )  and  horizontal  displacements  are 

denoted  u(t).  Thus,  the  solution  for  a  special 
case  of  Lamb’s  problem  may  be  written  as 
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and, 

S2  =  a  /  J3  ,y  =  {i  +  V3)12  /  2 , 
r  =  (j3/r)t, 

k2  =  [{a /  p)2  t 2  -  lj/[(a  /  ^)2  -  lj • 


The  k(&)  and  Yl(n,k)  are  complete  elliptic 
integrals  of  the  first  type  and  third  types, 
respectively. 


This  solution  proves  valid  for  Poisson’s  ratio  of 
v  =  0.25 .  A  solution  for  arbitrary  V  is  given  in 
[2] 


An  axisymmetric  model  of  a  semi-infinite 
domain  was  developed  in  COMSOL 
Multiphysics  to  solve  Lamb’s  problem.  For  both 
the  closed- form  and  finite  element  solution,  the 
forcing  function  is  described  as 


f(t)  =  for 

For  the  finite  element  analyses  conducted  in  this 
work  the  period  of  the  loading  is  10  jus, 

(T  =  10/^). 
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Figure  1.  Comparison  of  vertical  displacements  for 
a  three-dimensional  finite  element  solution  to  a 
closed-form  solution  for  point  loading  in  an  infinite 
domain. 


2.3  Comparison  with  Experimental  Data 

The  next  level  of  complexity  for  this  work  is  to 
represent  an  actual  geophysical  domain  over 
which  experimental  data  was  taken.  The  model 
developed  previously  in  COMSOL  Multiphysics 
was  modified  to  include  measured  pressure  wave 
speeds  (cp)  and  shear  wave  speeds  (c5)  and  to 
represent  the  actual  forcing  function  applied  the 
surface.  Over  a  depth  of  thirty  meters,  ten 
material  layers  that  represent  experimentally 
measured  values  of  cp  and  cs  values  were 
included  in  the  model.  The  pressure  wave 
velocity  increased  from  approximately  650  m/s 
to  2000  m/s  with  increasing  depth.  The  shear 
wave  velocities  increased  from  200  m/s  to  600 
m/s  over  this  range  of  depths. 

3.  Discussion 

3.1  Point  Source  in  an  Infinite  Media 

Error!  Reference  source  not  found,  shows  the 
comparison  of  vertical  displacements  over  a 
radial  line  from  the  center  of  the  sphere.  The 
results  show  the  transmission  of  the  wave  as  it 
moves  through  the  infinite  domain.  These  results 
show  the  ability  of  COMSOL  Multiphysics  to 
represent  accurately  the  wave  except  for  minor 
variations  at  the  peak  of  the  wave.  These 


differences  may  be  due  to  the  mesh  size  used  in 
this  model. 


Figure  2.  Comparison  of  radial  displacements  for  a 
three-dimensional  finite  element  solution  to  closed- 
form  solution  for  point  loading  in  an  infinite 
domain. 

Figure  2  compares  the  results  developed  for  an 
axisymmetric  model  of  an  infinite  domain.  These 
results  show  similar  accuracy  compared  to  the 
three-dimension  results  shown  in  Error! 
Reference  source  not  found..  Near  the  peak  of 
the  wave,  the  axisymmetric  model  shows  greater 
accuracy  than  the  three-dimensional  model.  This 
improved  accuracy  may  be  due  to  the  small  finite 
element  size  used  in  the  axisymmetric  models. 

This  work  also  used  a  plane  strain  formulation  to 
assess  the  ability  to  represent  a  point  source 
using  a  plane  strain  formulation.  Figure  3  clearly 
demonstrates  that  the  plane  strain  formulation 
provides  a  poor  approximation  to  a  point  source. 
The  plane  strain  formulation  transforms  a  point 
source  into  a  line  source.  Thus,  the  computed 
wave  form  changes  dramatically  between  the 
three-dimensional  and  plane  strain  solutions. 
These  results  indicate  that  an  axisymmetric 
formulation  provides  a  better  method  to  reduce 
the  computational  cost  of  this  problem  than  a 
plane  strain  formulation. 

3.2  Point  Source  in  a  Semi-Infinite  Media 

Figure  4  compares  results  developed  from  a 
finite  element  analysis  conducted  using 
COMSOL  Multiphysics  with  the  analytical 
solution  developed  in  Section  2.2.  These  results 


show  the  ability  of  the  finite  element  analysis  to 
represent  accurately  the  wave  propagation 
generated  from  a  volume  source  through  a  semi¬ 
infinite  domain.  The  finite  element  results  show 
the  lower  magnitude,  faster  moving  P  wave  that 
precedes  the  Rayleigh  wave.  The  finite  element 
analyses  accurately  calculate  the  arrival  time  and 
frequency  content  of  this  pulse.  However,  the 
finite  element  analyses  under  predict  the 
magnitude  of  the  velocity  by  approximately 
20%. 

3.3  Effects  of  Model  Details 

To  compare  with  experimental  data,  a  finite 
element  model  was  constructed  that  includes  the 
variation  of  the  soil  properties  as  described  in 
Section  2.3.  In  Figure  5,  the  effects  of  the 
variable  wave  speed  are  shown  compared  to  a 
model  that  has  homogenous  material  properties 
equal  to  the  first  layer  in  the  variable  model. 

The  homogeneous  model  shows  clearly  defined 
P  and  Rayleigh  waves  while  the  layered  model 
shows  the  wave  reflections  that  developed  due  to 
the  variations  in  material  properties  through  the 
depth  of  the  model.  The  shorter  arrival  time  for 
the  layered  model  represents  one  interesting 
feature  of  this  comparison:  the  increased  wave 
speed  below  the  surface  of  the  model  generates  a 
wave  that  reflects  off  a  subsurface  layer  and 
arrives  back  at  the  surface  in  less  time  than  the  P 
wave  traveling  in  the  first  layer  of  the  model. 

3.4  Comparison  with  Experimental  Data 

Figure  6  compares  the  radial  velocity  measured 
during  a  well-controlled  series  of  experiments 
with  results  from  the  finite  element  analyses.  The 
magnitude  of  the  velocities  has  been  normalized 
by  the  peak  magnitude  during  the  duration  of  the 
pulse.  This  normalization  was  done  to  provide  a 
comparison  using  the  two  primary  metrics  of 
interest  in  this  work:  arrival  times  and  frequency 
content.  Prediction  of  the  magnitude  of  the  wave 
is  a  secondary  consideration.  These  results  show 
a  reasonable  agreement  between  the  finite 
element  results  and  experimental  data  for  the 
primary  factors  of  interest. 


Figure  3.  Comparison  of  vertical  displacements  for 
three-dimensional  and  plane  strain  finite  element 
solution. 


Lamb's  Problem 


finite  element  solution  for  a  near-surface  point 
located  at  100  m  from  source  (T  =  10//s)* 


Figure  5.  Effect  of  varying  wave  speed  through 
layers  in  the  computational  model  (R=20  m). 


Figure  6.  Comparison  of  finite  element  results  with 
experimental  data  (R=20  m).  The  time  specified  on 
the  ordinate  is  relative  to  the  start  of  the 
experiment. 


4.  Summary  and  Conclusions 

The  results  presented  in  this  paper  demonstrate 
the  capabilities  of  a  COMSOL  Multiphysics  to 
solve  wave  propagation  problems  in  finite 
geophysical  domains.  The  analyses  were 
conducted  using  computational  hardware  that  is 
readily  available  on  the  desktop.  To  maximize 
the  use  of  desktop  hardware,  this  work  examines 
methods  for  reducing  the  problem  size.  For  the 
point  source  problems  considered  here, 
axisymmetric  modeling  provides  significantly 
more  accurate  results  compared  with  plane  strain 
modeling  techniques. 

This  work  initially  focuses  on  solving  problems 
in  a  uniform  infinite  space  or  infinite  half-space 
for  which  analytical  solutions  exist.  After 
showing  good  agreement  with  these  solutions, 
additional  complexity  is  added  to  represent 
available  experimental  data.  The  modeling 
methods  developed  in  this  work  again  show 
good  agreement  with  data  developed.  This  work 
shows  the  strong  effect  of  including  the  variation 
of  wave  speed  through  the  top  thirty  meters  of 
the  earth.  By  including  this  variation  in  material 
properties,  the  wave  arrival  time  and  frequency 
content  agree  with  experimental  data.  Thus,  this 
work  demonstrates  that  COMSOL  Multiphysics 
provides  a  useful  tool  for  predicting  wave 
propagation  using  commercially  available  finite 
element  software  with  desktop  computing 
hardware. 
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Overview 


•  Objective: 

Demonstrate  the  capability  of  COMSOL  Multiphysics  to  accurately  solve 
wave  propagation  problems  in  geophysics 

•  Motivation: 

Reduce  reliance  on  custom  software  and  supercomputers  by  obtaining 
solution  using  commercially  available  software  on  high-end  desktop 
computer 

•  Approach: 

-  Develop  closed-form  solutions  for 

•  Point  source  in  an  infinite  body 

•  Point  source  on  the  surface  of  a  semi-infinite  body  (Lamb’s  problem) 

-  Develop  using  solid  mechanics  module  w/  COMSOL 

•  Same  formulation  as  acoustics  module 

•  Three-dimensional 

•  Axisymmetric 

•  Plane  strain 

-  Comparison  w/  experimental  data 

•  Hammer  blow  on  surface 
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Closed  Form  Solution  -  Displacement 


Elastic  Wave  in  an  Infinite  Body 
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Point  Force  Solution  -  2D 


xlO 


2D  Model 
ND0F=  195,000 


•Axisym  metric 
•Plane  Strain 
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Comparison  of  3D  and  Axisymmetric 
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ACES 


Comparison  to  Analytical  Solution  -  Axisymmetric 
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ACES 


Comparison  of  3D  and  Plane  Strain 


Min:  -1.22 7e- 17 


Time  =8  Surface:  y-displacement  [m  ]  Mix:  7.663e-l5 

KlO"15 


Min:  -7.6756-15 


Advanced  Computational  &  Engineering  Services 


ACES 
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Surface  Wave  Problem 


2D  Model 
R  =  200  m 
Ndof  =  23,000 


Short  Duration  Loading 


Lamb's  Problem 
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ACES 
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Hammer  Data 
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Velocity  (m/s) 
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Summary 

•  Closed-form  solution  developed  for 

-  Elastic  wave  in  infinite  media 

-  Elastic  wave  in  semi-infinite  media 

•  Computation  models  developed  using  Solid 
Mechanics  Module 

-  Three-dimensional 

-  Two-dimensional 

•  Axisymmetric 

•  Plane  Strain  -  not  sufficiently  accurate  for  point  source 

•  Comparison  with  analytical  solutions  and 
experimental  data 

-  Agreement  with  arrival  time,  and  frequency  content 
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Conclusions 

•  COMSOL  Multiphysics  provides  a 
sufficient  level  of  accuracy  for  the 
problems  of  interest 

•  COMSOL  Multiphysics  provides  a 
commercially  available  tool  that  can  solve 
wave  propagation  problems  on  desktop 
computing  resources 
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